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Abstract 

We investigate static and dynamic properties of gray-scale image restora- 
tion (GSIR) by making use of the Q-Ising spin glass model, whose ladder 
symmetry allows to take in account the distance between two spins. We 
thus give an explicit expression of the Hamming distance between the orig- 
inal and restored images as a function of the hyper-parameters in the mean 
field limit. Finally, numerical simulations for real- world pictures are carried 
out to prove the efficiency of our model. 
PACS numbers : 02.50.-r, 05.20.-y, 05.50.-q 
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I. INTRODUCTION 



In the last decade, the problem of the image restoration (IR) has been successfully 
investigated by means of techniques borrowed from the field of statistical mechanics. 
Among them, it is certainly worthy mentioning the Maximum Posterior Marginal (MPM) 
and Maximum A Posteriori (MAP) estimations. From the statistical mechanical point of 
view, each recovered image within the MPM estimation can be regarded as the equilibrium 
state of ferromagnetic spin systems in the presence of random fields at finite temperature. 
In simple words, the reconstruction of a corrupted image is achieved by balancing the 
strength of a linear field, which carries the information of the degraded picture, and a 
ferromagnetic term which builds relatively large "one-color" clusters (below the transition 
temperature), thus suppressing the isolated pixels thought to be noise. From this point 
of view, the MAP estimation consists in the minimization of the same Hamiltonian at 
zero temperature (search for the ground state), with an appropriate scaling of the random 
field. The advantage of the MPM estimation over the MAP one has been pointed out 
by Marroquin et al [|l| and its performance has been investigated by several authors 
In this direction, Nishimori and Wong [Q, by unifying IR problem and error- 
correcting code theory under a single framework, found that the optimal recovering of an 
image is obtained at a finite temperature (known as Nishimori temperature in the field 
of statistical mechanics). Their results, however, were restricted to the usual binary spin 
models (Ising), i.e. black or white images in IR jargon, and many questions about the 
properties of the gray-scale image restoration (GSIR) processes still remain open. A first 
attempt to generalize |4j to gray-level pictures has been carried out by the authors || 
by mapping the set of the pixels onto Q-state (chiral) Potts spins in the presence of the 
random fields. In that case, the symmetry of the Potts Hamiltonian {hyper-tetrahedron) 
reduces the problem to a 2-stoXe-like system, where only one bit turns out to be right, 
and all the others are equivalently wrong without any regard of the whole gray-level scale. 
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Whereas this turns out to be an efficient method in the presence of white noise (each spin 
is flipped to any of the Q values with equal probability), things may be different from 
a transmission channel affected by Gaussian noise (the spin-flip probability distribution 
is a Gaussian). We thus investigate the performance (both static and dynamic) of the 
gray-scaled image restoration using the Q-Ising model 0, whose ladder symmetry takes 
in account the distance between the spin values and will allow us to say, for instance, 
that Q = 3 is better than Q = 5 if the right pixel corresponds to Q — 2. The analytical 
expressions are obtained in so-called mean field limit, where each spin interacts with 
all the others. The efficiency of our model is checked by Monte Carlo simulations and 
iterative algorithm by using mean-field approximation. This paper is organized as follow. 
In the next section, we introduce the infinite range Q-Ising model and in Sec. |j A| we give 
an analytical expression of the Hamming distance in the mean field limit. In Sec. |1 TB| 
we derive the dynamical equations with respect to the macroscopic quantities, namely, 
the magnetization and the Hamming distance in terms of microscopic master equation. 



In Sec. |TJ, in order to test the usefulness of the Q-Ising model for the GSIR, we carry 
out Monte Carlo simulations for real-world pictures with Q = 8 gray-scale levels. In 



Sec. [TV], we show an iterative algorithm based on the mean-field approximation, whose 
convergence is much faster than that of the Monte Carlo simulations. The last section is 
left for summary and discussions. 



II. THE INFINITE RANGE Q-ISING SPIN GLASS MODEL 



A Q-gray scale levels image is nothing else that a set of pixels {£} on a grid, whose 
values can be coded at each node as an integer variable G {1, 2, • • -, Q}. Without loss 
of generality, let our image be generated by the following prior distribution 



P s ({£}) = — exp 



where Z s is the usual normalization constant that is given by 



Z s = tr^exp 



ij 



(2) 



In the spirit of statistical mechanics, we want to regard this picture as a snapshot of a 
spin system described by the Hamiltonian H s = (1/2) J2ij(£,i — £7) 2 at a specific temper- 
ature T s = /S^ 1 - Sending our image through a noisy channel will cause the flipping of 
some pixels to different values. For this degrading process, we assume that each pixel & 
changes its state to 7$ independently. Then, the degraded pixel r, is given by the following 
conditional probability 



1 



2vrr 



exp 



2r 2 



(3) 



This means that after the transmission, the receiver observes Tj which was violated from 
scaled original image r £j with a standard deviation r. This kind of damaging process 
is referred to as "Gaussian channel (GC)". Due to the independence of noisy process on 
each pixel, a sequence of original pixel {£} is corrupted by the GC as 



i 

1 

exp 



2ttt 



2r 2 



Efr-T0fc) 5 



(4) 



In the context of Bayesian approach, the probability that the estimate of the source 
sequence is {a} provided that the observed noisy data is {r}, reads 

P({r}\W})P({a}) 



P(W}\{r}) 



tr W} P({r}\{a})P({a}) 



exp 




hT.M-nf-^ 






tr M exp 


-hEMi-Ti) 2 - 


(A*/2) - 



_ exp (-Weff) 



(5) 



where we defined the effective Hamiltonian H e s and the normalization constant as 

(3d 



Ties = hJ2(ai - Tif + — J2( a i - G 3 



(6) 
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and 



Z d = tr {(T} exp [-Ties] 



(7) 



respectively. The parameters h and f3 d appearing in the Hamiltonian 7i c ^ [Eq. (^|)] are 
referred to as "hyper-parameters" and we can not mention about the true values of them 
beforehand. This conditional probability -P({c"}|{t}) is called "posterior probability" and 
is constructed in terms of a likelihood -P({t}|{c"}) and a prior probability P({a}) as we 
saw in Eq. (^). -P({t}|{<t}) and P({a}) are given by 

exp [-hJ2i{n - Pi) 2 } ,„s 



^(MIW) 



tr {T} exp [-hT,i(n - o { 



and 



p(W}) 



exp 


-W2)Eij(o-i-o-i) 2 " 




tr M exp 


-W2)E^(^-^) 2 



(9) 



The prior probability reflects our assumption on the original image that the picture should 
be locally smooth. As shortly mentioned in the introduction, the MAP estimation consists 
in maximizing the above posterior probability P({a}|{r}), that is finding the ground state 
{a} of the effective Hamiltonian 7i e g and regarding it as an estimate of true pixels. 

On the other hand, in the context of the MPM estimation, we first consider the 
following marginal distribution; 



Pfali*}) = E P(W}\{r}) 

and then we calculate the local magnetization which is given by 

Q 

<Tj = l 



(10) 





exp 


-hYUfa - a i) 2 - (/V 2 ) - 




tr {(T }exp 


-h Eiin - °i) 2 - (f3 d /2) Zijfa - o^) 2 " 





tr^c^exp 


-hJ2i(Ti - (Ti 


) 2 -((3 d /2) ZiM-°i) 2 '_ 


tr M exp 


-hYAn - ViY 


1 ~ (fid/2) Yij(°i - a j) 2 


tr{ CT }(Tiexp 


[— Hes] 




tr {(T} exp [ 


—T~tes] 
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Using the above expectation value, we regard the estimate of the original pixel & as 
^{( a i) Pd,h) where function Q is represented by a sum of step functions 0(x) (Q(x) = 1 
for x > and 0(x) = for x < 0); 



fc=i 



{*i)f)i,h. 



2k - 1 



<^i>/»d,h ~ 



2fc + 1 



(12) 



2 j y IHd) 2 , 

The natural quantity measuring the quality of our restoration process, viz. the distance 
between the original and the recovered image, is the Hamming distance (square error) 



(13) 



whose value depends upon the hyper-parameters, h, fid appearing in the effective Hamil- 
tonian 7i c ff. At this stage, it is important to bear in mind that the MAP estimate is 
recovered as the limit >oo (keeping their ratio constant H = h/fta) in Eq. (|TT|). En- 
couraged by the results in Q and |5[], we expect that more data fed through the noisy 
channel improves the quality of the restored image, since the receiver will have more infor- 
mation about the original image. Therefore, in addition to the transmission of the single 
bit £j, we send also the pairwise product Then, each product is also corrupted 
independently by the following GC 



2nJ 



exp 



2 J 2 



(Jij — Jo£t£j 



(14) 



namely, the degraded version of the product deviated from the scaled original data 
JoCiCj with width J. For this degrading process, we modify a likelihood -P({t}|{o"}) in 
Eq. (D as 



P({t},V}\{*}) 



exp 



(A// 2 ) J2ij(Jij ~ cria 3 ) 2 -hYM - Oi 



(15) 



with 



z l = tr {r},{j>exp 



h^2(Ti - <7if 



(16) 
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where we introduced another hyper-parameter f3j. Using the same way as Eq. 
rewrite the posterior probability as 

exp [-H e s] 



P({°}\{t},{J}) 



tr {(T} exp [-H c $] 



we 



(17) 



with the following effective Hamiltonian 
T^eff = Jii — OiOj) 2 — h 



2 ^ 



Ti — (Ti 



Given the degraded version of data, namely, {r} and {J}, arbitrary macroscopic phys- 
ical quantity f({u}, {r}, {J}) is calculated in terms of the average over the posterior 
distribution P({cx}|{r}, {J}) as 



</(W, to, {J}))0 d ,h = tr W /(M, {r}, { J})P(M|{r}, {J}) 

_tr w /(W,{r},{J})e-^ 



(19) 



tr{ CT }e Wrff 

As the quantity (/({o"}, {r}, { J}))/3 d ,h depends on the observed data {r}, {J}, we should 
average them out by the distribution 

P({r},{J}\m) 

(20) 

Therefore, after tracing the original image {£} out, the averaged macroscopic quantity is 
given by 



exp 


i 


EdJij - - (i/2r 2 ) Eiin - &) 2 - f Zijfa - & 2 _ 




tr M,{J},{«>exp 


^a( J a ~ U£i) 2 - (V2r 2 ) J2 t (n - e t ) 2 - f E»;(& - 


- Q\ 



K/(M> to, { J }))p d ,h}{r},{J},{0 

"tr W /(M, to, {J}) e -*- 

tr W,{^},{«} 



p{{t},{J}\{Z}) 



(21) 



tr w e-^eff 

Using this definition, the performance of image restoration is measured by the following 
averaged Hamming distance between the original image and the restored one, that is, 
to((<7i)fi d ,h) as 

-tr {a} (l/2N)j:Mi-n((°i)(S d ,h)) 2 e 



DK = ^{r},{J},{0 



tr{ CT }e HcS 
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P({r}AJ}\m (22) 



In the next two subsection, we investigate the performance of image restoration in terms 
of this Hamming distance _D#. We focus our analysis not only on the static properties 
but also the dynamic properties of image restoration. 



A. Static properties 

In this subsection, we consider the static properties of image restoration. First of all, 
we should investigate the properties of original image, that is to say, the properties of the 
ferro magnetic Q-Ising model. However, it is quite hard to calculate the partition function 
or the other physical quantities for our spin system defined on two dimensional square 
lattice analytically. Therefore, in this paper, we investigate the infinite range version of 
our model system and calculate the macroscopic physical quantities analytically. Then, 
the infinite range version of the prior distribution leads to 



F.({£}) = ^exp 



2N 



(23) 



where we should notice that the argument of the exponential should be divided by N in 
order to take a proper thermodynamic limit. For this rather artificial model, we easily 
obtain the magnetization at some temperature T s (= ft' 1 ) as follows. 



N 

tr ? ^e 2mo/3s€ ~ /3 ^ 2 



(24) 



We should be in mind that for the infinite range Q-Ising model, the properties of the 
macroscopic quantities of the system are completely determined by m . In FIG. [1], we 
plot the magnetization mo as a function of source temperature T s for Q = 3 and Q = 4. 
We see that for the Q = 3 case the three states m = 1,2 and 3 are degenerated at 
T s = 0, while at finite temperature the middle state mo = 2 becomes globally stable and 
mo = 1,3 are degenerated locally stable. At high temperature regime T s — > oo, each spin 
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takes all the values with same probability 1/3 and thus the corresponding magnetization 
is m = (1 + 2 + 3) /3 = 2. The transition between the ferro- magnetic phase and the para- 
magnetic phase occurs at T c ~ 1.0. In the same way as the case of Q = 3, for Q = 4, the 
four states mo = 1, 2, 3 and 4 are degenerated at T s = 0, and the middle two states itlq = 2 
and 3 become globally stable for T s > (m = 1,4 are degenerated locally stable states). 
The para-magnetic state is specified by the magnetization m-o = (1 + 2 + 3 + 4)/4 = 2.5 
and the ferro-para transition occurs at T c ~1.78. For this original image, in order to 
investigate the average performance of the MPM estimation, we should calculate Dh in 
terms of statistical mechanics of the spin system {a} with quenched disorder {r}, {J} 
and {£}. For this purpose, we calculate the averaged free energy of the system described 
by TC e ff [Eq. flilf) ] with assistant of the replica method; 

[logZj Mi{J}j{€} - lim (25) 

with 

Z = tr w exp(-/?Kfr), (26) 

which consist in replacing the quenched average of a single system with an annealed 
average of n replicated systems (letting n — > at the end). Assuming a replica symmetric 
ansatz and by using the saddle point method, the order parameters are given by the 
following coupled equations; 

/oo 
VxB(x,0 (27) 
-oo 

/oo 
VxB(x,£) (28) 
-oo 

/■oo 

Q = M)^i)^h]M,{JUr} = tr c e(£) / Vx [B(x,0] 2 (29) 

J — oo 

/oo 
VxC(x,0- (30) 
-oo 

where we should remember that the brackets (• • -)p d .h and [■ • ']{r},{j}.{§} are defined by 
Eq. ( |i~9l) and Eq. fl2l|) , respectively. In the above expressions, Vx = (dx / y/2n)e~ x2 / 2 is 
the usual Gaussian measure and we defined 
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B(x,0 



tr^crexp [Ua — Va 2 } 



(31) 



C(x,0 



tr CT exp [Ua — Va 2 } 
tr a a 2 exp [Ua - Va 2 } 



(32) 



2(0 



tr CT exp [Ua — Va 2 } 



(33) 



with 



U/2 



{r h + (3jJ t)C + mf5 d + x^r 2 h 2 + (3 2 J 2 q 



(34) 



V 



h + (3 d + /3jw. 



(35) 



Using the same way as the derivation of the order parameters, the averaged Hamming 
distance Eq. (^) is calculated and reads 




case. To keep things easy, we first assume that there is no glassy term in our decoding 
process, i.e. (3j = 0. In FIG. 2, we plot the Hamming distance as function of the decoding 
temperature for Q = 3 and T s = 0.75. The minimum is reached at T d = 0.75. The same 
for Q = 4 in FIG. 0. It can be shown numerically, at least for Q = 3 and Q — 4, that, 
given the original image at temperature T s , just below the transition temperature, the 
optimal decoding temperature Tj opt ' ) = T s . The same relation turns out to be satisfied 
for black or white IR [|J (which corresponds to Q = 2), differently from GSIR by the 
Potts model 0. In order to compare the performance of the MPM estimate with that 
of the MAP one, we first investigate the scaled field H = h//3d dependence of the MAP 
estimate. The MAP estimate is obtained by controlling the temperature as T d — > with 
keeping the scaled field H constant. Therefore, the Hamming distance for the MAP 
estimate should depends on H . In FIG. 4 (a), we plot the Hamming distance of the MAP 
estimate as a function of H. In this figure, we set Q = 3, T s = 0.75 and r = tq = 1.0. We 
see that the Hamming distance takes its minimum at H opt = t /2t 2 (3 s = 0.375, namely, 





(36) 



It is straightforward to check that the above equations coincide with those in for Q = 2 
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Du(T d = 0, H) > Dn(T d = 0, t /2t 2 /? s ). This optimal value of the scaled field H = h//3 d 
is obtained when we set -P({t}|{£}) = -P({ r }|{ cr }) an d = PdiW}), that is to say, 

Pd = Ps and h = r /2r 2 . In FIG. 4 (b), we increase the temperature T d with H = H opt 
and plot the Hamming distance Dn{T d , H opt ) as a function of T d . This figure shows that 
Du(T d ,H opt ) takes its maximum at T d = T s = 0.75. Therefore, we conclude that the 
MPM estimate achieves the lowest Hamming distance which can not be obtained by the 
MAP estimation. In FIG. 4 (b), we plot the D}i(T d ,H) for several values of H. From 
this figure, we see that as long as we choose H so as to satisfy H > H opt = t /2t 2 (3 s , the 
minimum value of the Hamming distance does not change. 

In the limit of T d — > oo, each pixel takes <7j = 1,2,3 with same probability 1/3, and 
the local magnetization leads to (cr^) = (1 + 2 + 3)/3 = 2 for all pixels i. As the result, 
the Hamming distance in the high temperature limit takes 

D u (T d -> oo) = ^~ 1V 5 i nc nc2 0.1726. (37) 

This asymptotic behavior is checked in FIG. 4 (b). We now switch on the product inter- 
action, that is, (3 j 7^ setting T d and H at their optimal values. As clearly shown in FIG. 
[|, the performance of the restoration is dramatically improved. In this figure, the point 
(3 j = corresponds to the minimum of the Hamming distance in FIG. [|. 

B. Dynamics 

An important and interesting problem is to determine the basin of attraction of the 
Hamming distance Du(t). In fact, because of the presence of locally stable states, the 
final state of the decoding process is strongly dependent upon the initial condition of the 
dynamics. In addition, as the number of local minima increases with the the number of 
the gray-scale levels, it becomes crucial to choose the initial state appropriately. However, 
as it is well known, it is difficult to treat the dynamics of spin system explicitly in finite 
dimension, especially, dynamics in two dimension which is the case of image restoration. 
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In the previous section, we introduced the infinite range model and solved it analytically. 
Using this model, we derived the properties of image restoration and as we see in the next 
section, the results do not contradict qualitatively with the properties in two dimension. 
With this fact in mind, we also use the infinite range model to investigate the dynamical 
properties of image restoration. In the equilibrium limit t — > oo, without glassy term, 
namely, (3j — 0, the properties of image restoration in the infinite range model are com- 
pletely written by magnetization m. Therefore, we assume that the dynamics of image 
restoration is also expressed by the time evolution of the magnetization m(t). Therefore, 
we derive the differential equations with respect to the macroscopic variables, namely, 
m(t) and Dn(t), from the microscopic master equation. For the sake of simplicity, we 
restrict ourselves to the case without the glassy term. The master equation of our system 
leads to 

dptti*}) 



dt 



N Q 

E £ [ w ( a k' Vk)pt(W}') - w(a k -> oy)p t ({o-}) 



(3* 



fc=l a. i=l 

k 



with the following transition probability 



-(h+f3 d )a 2 , +2(m/3 d +hT k )cr , 

Q k « 



-(h+p d )a 2 ,+2(mt3 d +hT k )o,' 

k 



(39) 



where we defined {a} = (<7i, ■ ■ -, a k , • • •, <?n) and {a}' = (o"i, ■ • •, a k >, ■ ■ •, a^). By intro- 



ducing the probability distribution of the macroscopic magnetization m ,viz. 

'Pt{m) = ^Pt{{o})b[m - m({a})) 

a 

and after some algebra we obtain the following exact differential equation 

Q /-OO 

-m + £Q(£) / DxB(x,£) 



(40) 



dm 
~dt 



(41) 



The derivation of the above differential equation is reported in Appendix A. The time 
evolution of the Hamming distance -Dh(^) is obtained by substituting the time dependence 
of the magnetization m(t) into Dn(m). In FIG. [], we plot the time evolutions of the Q = 3 
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Hamming distance for Td = T s (a) and Td 7^ T s (b). From these figures, we see that if 
we choose the hyper-parameter Td so as to satisfy the relationship Td = T s , the Hamming 
distance converges to its optimal value for any initial condition. On the other hand, 
for Td 7^ T s , there exists a threshold of the initial value of the Hamming distance 
beyond which the flow Du(t) does not converges to its optimal value. As the dynamical 
equation (with respect to m) is exactly the same as the Time Dependent Ginsburg-Landau 
(TDGL) equation, that is, dm/dt = —dfus/dm, the nature of the dynamics is intuitively 
understood as a steepest descent to a local minimum of the free energy. In fact, from 
FIG. (a),(b) and (c), we see that for Td < 0.35, there exist local minima. Therefore if 
we fail to choose the initial condition appropriately, the Hamming distance converges to 
the non-optimal values. For practical situations, the corrupted image corresponds to our 
initial state. For the case Td 7^ T s , we calculate the Hamming distance D$ between the 
original image and the noised one, which reads 



nW = _L Vf 



n - & 



2Z{(3 S 



(42) 



In particular, for r = r = 1.0, this leads to D# = 0.5. From FIG. 5 (b), we see that 
if we choose the corrupted image as an initial state, we destroy the observed corrupted 
image and the result is even worse. 

The asymptotic expressions of the magnetization and the Hamming distance in the 
limit of t — > 00 lead to 



m = + [m(t = 0) — m*] e 'o . 



(43) 
(44) 



where m* is a solution of the saddle point equations Eqs. 
previous section. The relaxation time to is given as 



- = 1 + 2/3 d 

to ■' - x 



a- 



(a) 2 



with 



D with pj = in the 



(45) 
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\ e 2(m t l3 d +Thx+T h^)u~(h+l3 d )a 2 



y^V p 2(m»/3 d +T/ia;+To/i5)o--(h+/3 ( i)o- 2 

^(7 = 1 ^ 



and i)n reads 



(46) 



D H = A(3 d [m{t = 0) - m,] / Dx[£ - Q(w)] 



x 



fc=i 



5 



2k 



a 



5 



2k 



a 



(47) 



We plot the inverse relaxation time 1/to as a function of for the case of Q = 3, T s = 
0.75, tq = t = 1.0 in FIG. 7. In this figure, we also plot the inverse relaxation time 
for several values of the scaled field H. We see that the inverse relaxation time 1/to 
takes its minimum at a finite temperature However, the inverse relaxation time 1/to 
never reaches to zero, and the relaxation to the equilibrium state is exponential for all 
temperature T d regions. 



III. MONTE CARLO SIMULATIONS 

So far, we worked under the assumption that all the pixels lay on an infinite dimen- 
sional grid, an approximation which enabled us to derive exact analytical formulas. In 
order to test the efficiency of the Q-Ising model on the more realistic case of two dimen- 
sional picture, in this section, we carry out Monte Carlo simulations at finite temperature 
on a real-world image with Q = 8 gray-scale level [FIG. |9] (a)] corrupted by a Gaussian 
noise with r = 1.2 [FIG. || (b)]. Here the interaction in effective Hamiltonian is now 
restricted to the nearest neighbors spins on two dimensional square lattice. As before, we 
first study the Hamming distance without the glassy term. The resulting curves averaged 
over twenty Monte Carlo runs are shown in FIG. (^j) (a) for three different values of the 
ratio H = h/fta- The plots reflect indeed the mean field behavior of FIG. |^ (d) and [3] 
(b). The corresponding restored picture at optimal values is shown in FIG. |9| (c). It is 
evident that the ferromagnetic term succeeds in eliminating the noised pixel, i.e. isolated 
ones, but at the same time it also smoothes out the small true details of the original 
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picture. For this reason, by keeping fixed T d and H opU we switch on the glassy term, 
namely, (3j 7^ 0. Here we choose a slightly more general expression for Jy = (£j — £j) 2 
and therefore the extra term reads ([J^ — (a.; — o\,') 2 ] 2 /iV). The curves in FIG. |8] (b) for 
twenty Monte Carlo runs, show an improvement of the recovered image. The restored 
image at the minimum /3j° pt ' ) is on the lower right corner of FIG. |S|. It is evident that 
the extra term preserves many of the small details of the original image, for example the 
white edge of the roof, which was blurred for the case of (3j = 0. 



IV. ITERATIVE ALGORITHM (MEAN-FIELD) 

The restoration by means of Monte Carlo methods is the result of a statistical process 
which might take long time even for powerful computer as the size of the picture increases. 
Therefore we apply an iterative algorithm, proposed in |||9[], to our model. Using the 
mean-field approximation with periodic boundary conditions for two dimensional square 
lattice of size L\ x L 2 , the recursion relations with respect to the local magnetization at 
a site (i, j), namely, rriij lead to 

- , (48) 

tr CT e^>) 

4f(a) = { J[mg+i + ™S-i + mfl 1;j + mfl 1;j ] + 2hr l3 } a - (2J + h)a\ (49) 

m it j + L 2 = rriij, rn i+Lli j = m^. (50) 

where tr cr (- ■ ■) means the sum with respect to the gray scale levels, namely, a — 1, 2, - • -, Q. 
The details of the derivation by using a variational principle are reported in Appendix B. 

Then, we obtain the estimate of the pixel namely, f^m^) by solving the above 
non-linear maps until appropriate error tolerance is satisfied. In order to investigate its 
performance, we introduce the following three measures 
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^H = ^VEE[fiK,)-^] 2 (si) 

1 z J = l J = l 

4 1} = E f>M - r,,] 2 (52) 



1 z «=i j=i 



where , and are distances between the original image £y and the restored 
one ^(mij), the corrupted image and the restored one, the corrupted image and the 
original one, respectively. We choose Q = 8, L\ = L2 = 128 and h = 1.0 and solve the 
recursion relations (]48l),(]49f) and (1501) until the error 



^-^VEEkr-^l (54) 

^1^2 i= i j= i 

becomes smaller than 10~ 5 . We list the results in FIG. 8 ("girl") and FIG. 9 ("chair"). The 
original images are degraded by the Gaussian noise with a standard deviation r = 2.2. 
This standard deviation gives the Hamming distance between the original image and 
degraded one -Dh 0.205. Obviously, the picture "chair" contains much more edges 
than the picture "girl" . Therefore, one of the our aims of this demonstration is to check 
to what extent our model can detect the edge parts in the real-world picture. In FIG. 
10, we plot the Hamming distance Du [(a)] and Djj [(b)]. We choose the degraded 
image as an initial set of the pixels and investigated the J-dependence of the Hamming 
distance. We see that the performance of the algorithm for the "girl" picture is much 
better than that of the "chair" picture. This is because the smoothness term in the 
effective Hamiltonian [Eq. ([])] is quadratic and it is hard to detect the edges in the "chair" 
picture. For both pictures, the optimal performance is achieved around the parameter 
Td = 1/J = 1/1.8 ~ 0.56. This value is not so different from the parameter which was 
obtained in Monte Carlo simulations [see FIG. 6 (a)]. Of course, from a practical point of 
view, it is possible to stop the Monte Carlo simulation and not to wait the convergence 
to the equilibrium state precisely. Then, we may regard the snap-shot as the restored 
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image if the performance is not so bad. However, in the mean-field approximation we 
constructed here, the convergence of the iteration is guaranteed. 

V. SUMMARY AND DISCUSSIONS 

In this paper we investigated the efficiency of the Q-Ising model for image restoration 
problem, when the original image is affected by Gaussian noise. By introducing the 
infinite range model, we gave an analytical expression for the Hamming distance, which 
is shown to reach its minimum at some finite temperature. We found that the optimal 
temperature for the GSIR using the Q-Ising model coincides with the source temperature 
in contrast to the chiral Potts case ||. We also found that as in the Ising and Potts 
spin cases, the presence of a parity-check- like term greatly increasing the performance 
of the GSIR process. Although for practical restorations of images, one wouldn't like to 
smooth out two points far away, the mean field results provide a remarkable eye-guide for 
a short-range version of the effective Hamiltonian as confirmed by Monte Carlo simulation 
on two dimensional pictures. From a dynamical point of view, we also obtained the time 
evolution of the Hamming distance analytically and found the critical initial Hamming 
distance beyond which the flow does not converges to its optimal value. We show the 
dynamical equation we obtained is exactly same as the TDGL equation. Therefore, the 
destination of the dynamics is one of the local minima of the free energy, and if we fail to 
select the initial condition, the dynamics converges to the local minimum which does not 
give the minimum of the Hamming distance. 

Recently, Skantzos and Coolen flC| , reported the synchronous dynamics of 1-D and 
infinite range version of the random field Ising model. They found that the dynamics 
has much more rich behavior than the sequential (Glauber) dynamics. Therefore, for our 
present model system, there is a possibility that if we consider the synchronous dynamics 
instead of the sequential one, the behavior of the dynamics may be different from the 
results we obtained here. 
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Using the mean-field approximation, we also constructed the iterative algorithm which 
converges faster than the Monte Carlo simulation. We derived it from a variational 
principle of the free energy and demonstrate it for two types of the real-world pictures. 
From those results, we concluded that we need some extra term which detect the edges 
if the picture has a lot of edges. We suppose that the glassy term we introduced in the 
infinite range model may play this role. This will be achieved by means of the TAP like 
mean-field approximation. 

The authors acknowledge Profs. Hidetoshi Nishimori, Kazuyuki Tanaka, and Desire 
Bolle' for fruitful discussions and useful comments. We also thank Dr. Masato Okada for 
useful discussions about the dynamics of disordered systems. 

One of the authors (J. I.) was partially supported by the Ministry of Education, 
Science, Sports and Culture, Grant-in-Aid for Encouragement of Young Scientists, No. 
11740225, 1999-2000. 

APPENDIX A: DERIVATION OF THE FLOW OF MAGNETIZATION 

In this Appendix, we derive the differential equation with respect to macroscopic order 
parameter m from microscopic master equation for the infinite range version of the Q- 
Ising model. For simplicity, we consider the case of no-parity check term j3j = 0. For the 
Q-Ising model, the effective Hamiltonian is given as 

= ^ - ° 3 ) 2 + hJ2(^ ~ nf = H(a k ). (Al) 

ij i 

Therefore, the energy difference due to the local spin change a k — > a k > , namely, AE = 
7~t(a k >) — H(<j k ) is calculated in terms of the above Hamiltonians H(a k ) and H(cr k ') as 
follows. 

AE — (h + (3 d )(a 2 k , - o\) - 2(3 d m{a k > ~ " M°k' - °k)r k , (A2) 

where we used the expression of the magnetization 
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N A.. 



(A3) 



Then, the transition probability w(a k — > ay) is given by 



-AE 



w(a k ^a k >) = —q 



k 



AE 

(h+/3 d )a 2 , +2(m/3 d +hT k )a, 



srQ ~(h+/3 d )a 2 +2(ml3 d +hT k )a k , 

E V =l e 



For this transition probability, the master equation leads to 

d n Q 

^Pt(M) = E E [ w ( a k' W') - ^(^ -" <v)Pt(M) 



fc=l <T ,=1 

fe 



Here we introduce the following macroscopic probability 

M 

and consider the derivative of Vt(m) with respect to i, that is, 



(A4) 



(A5) 



(A6) 



V t (m) = ^ 6[m - m{{a})\ 



dt 



M 



dt 



( N Q 

E E E l w (°k'^k)pt(W}) - w(a k ^a k ,)p t ({a})} \ 5[m - m({a})} 
W U =1 V =1 

J2 w ( a k^ k ')Pt(W}) U m - m({a}) + — (a k - ay) - 5[m - m({a})\ 
Mfc=i v =i L 

Q ( N Q 1 

a~ \Y,Y1 £ w ( a k^k')Pt(W})^[rn-m({a})} — (a k -a k ') 

° [W}k=la k ,=l 

d { N Q 



-(/t+/3 d )<r 2 , +2(m/9 d +nfe/i)<7 / 

Q k K 

Q -(/i+/3d)^, +2(m(3 d +T k h)a k 



fc=l ^ (7, / =1 



Q 

- V 

iV 

fX / =1 



k =1 

-(/i+/3 d )<r 2 , +2(mf5 d +T k h)a , 
v^Q -(h+p d )o 2 ,+2(m(3 d +T k h)o k , 



-{h+/3 d )<r 2 +2(mP d +T k h)a , 

-(/i+/3 d )a2 +2(m/3 d +r fe h)<7 fe 



>£[m — ?n({cr})] 
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{a} y 

x 8[m — m({a})]. 



N 



M v^Q -(/i+/3 d )(7 2 / +2(m/3 d +T fe h)CT 
fc ^ v^Q -(h+/3 d ) CT 2 ,+2(m/3 d +T fc /iW , 

fc=i fc=i e 



(A7) 



Here we should notice that 



-(h+/3 d )cr 2 , +2{ml3 d +Th)a k , 



iV V-Q Q -(ft+/3 d )<7 2 ,+2(m/3 d +r/ l W , 

V-i 



2(A) 



"Da; 



-(/i+/3 d )o- 2 +2(m/3 d +hTX+hro?)cr 



£ ( 3_ i e _(/ l+/ g d ) CT 2 +2 ( m /3 d+ ?j ra ; + / l713 ^) <T 



(A* 



should be hold due to the self- averaging properties in the thermo-dynamical limit iV — > oo. 



Substituting this expression into Eq. (|A7|) , we obtain 



d d d 

V t {m) = — m^2p t ({a})5[m - m({a})} - — ^p t ({a})5[m - m({cr})] 

M {<?"} 

2 



x 



^«_ ie 2mo/3 s 5-/3 s C 2 



rriPt(m) 



Vx 



ae -{h+l3 d )o 2 +2{ml3 d +hTX+hT £,)<j 
Q-{h+P d )a 2 +2(mj3 d +hTX+hT S,)u 



_8_ 

dm 



d ,E?=ie 



Q „2m /3 s £-/3 s £ 2 



x, Tv-Q -{h+j3 d )o 2 +2{m0 d +hTX+hT O £,)a 

Vx ^ a=1 

oo I Y^ = iQ~^ h+Pd ^ 2+2{mf3d+hTX+hT0 ^ a 



-^—{v t (m) \ in 
am I 



(•OO 

x / Vx 

-oo 



Ect-i cre _ ^ +/3d ) <T2+2 ( m/3d+ ' lTX+/lro?)CT 



(A9) 



Multiplying m and substituting Vt(m) = S[m — m(t)} to the left hand side of the above 
Eq. (ggfl , we obtain 

/ mdm— 5[m — m{t)} = — ( m dm 5[m — m{t)] = — — . (A10) 

J —OO (it (It J —OO (it 

Using the same way as the left hand side of the Eq. (]A9|) , the right-hand side of Eq. (|A9|) 
leads to 



mdm— — 
dm 



5[m — m{t)] ( m 
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x, Ty-Q -{h+f3 d )a 2 +2(m[3 d +hTX+hT Z)a 

Vx ^ a=1 

oo [ Y^_ l e-( h+ P^ a2+2 ( m i 3d+hTX+hT0 & cr 



dm S[m — m{t)] < m — 



Vx 



Z{(3 S ) 

J2 Q _ 1 ae~( h+l3 ^ a2+2 ( ml3d+hTX+hT0 ^' 7 

Y^ =ie -(h+t3 d )a 2 +2(m(3 d +hTX+hT Z)a 



—m + 



^Q_ ie 2m /3^-/3 s ? 2 



Vx 



Y^_ l(Te -( h +Pdh 2 +2(mf3 d +hTX+hTo£,)<T 



X^Q p-(h+(5 d )a 2 +2(m0 d +hTX+hT O i)c 



(All) 



From Eqs. (|A10 ) and (|A11|) , we obtain the final form of the dynamical equation with 
respect to magnetization m as 

E?-ie 2m ° /3s5 ~ /3s?2 f°° r^Q_ i0 - e -(ft+/3 d )fT 2 +2(m/3 [i +/ l TX+/iroO^" 



dm 
~dt 



-m + 



Vx 



(A12) 



J2^ =ie -(h+l3 d )a 2 +2(ml3 d +hTX+hT ^a 

We easily see that the above equation is exactly same as the Time Dependent Ginsburg- 
Landau (TDGL) equation which is derived from a steepest descent of the replica symmet- 
ric free energy, that is, —dfus/dm = dm/dt. We should also notice that in the limit of 
t — > oo and dm/dt = 0, Eq. ( |A12| ) corresponds to the saddle point equation with respect 
to m which was calculated by equilibrium statistical mechanics in Sec. IIA. 



APPENDIX B: VARIATIONAL PRINCIPLE FOR THE Q-ISING MODEL 

In Sec. VI, we introduced the recursion relations which determine the estimate of 
the original image in terms of mean-field approximation f||J. In this appendix, we show 
that these recursion relations Eqs. (f48l) , fl49"l) and fl50|) can be derived from a variational 
principle. 

We consider the following optimization problem; 

min p {£(p)-TS(p)} (Bl) 

^)=EW(W)/>(W) (B2) 
M 

^)^-Ep(W)MW), (B3) 

W 
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where Hamiltonian Ti is defined on the two dimensional square lattice of size L\ x L 2 
(L 1 = L 2 = N) as 

7~L{{(y}) = - o"ij+i) 2 + - °i,j+i) 2 + - a i+i,j) 2 + (cy - CTi-l,i) 2 } 

y 

+ ^E(^-^) 2 > (B4) 

y 

and £ and 5 correspond to the energy and entropy of the system, respectively. Then, we 
use the mean-field approximation, that is, 

p(W) * \\l>.Mo)- (B5) 

ij 

We should notice that for each pixel the following normalization condition should 

hold 

Q 

E Pufa) = 1- (B6) 

<7y=l 

Using the Lagrange multiplier Ay, we take into account the above normalization condition 
with respect to the marginal distribution, and maximize the following functional 

T = £({a}) - TS{{a}) + E Ay ( E Pn{°*) ~ 1 ] • ( B 7) 
The energy £ and the entropy S of the system can be written explicitly as 
£ = E""E"'E"' E W({<t})pi 2 (<7i 2 )- • -Py(o-y)- • -Pki{crki)- ■ -Pn-in(<tn-in), ( b 8) 

CT 12 Cjj ""fci (TN-1N 

E n^K)E ln ^K)- ( B Q) 

""12 T13 (Tjj 0"jV-ljV ij ij 

The derivative of the third term of Eq. ( [B7D with respect to /0y(cy) leads to X^ J= i Ay, 
therefore, we have [dT jdpij) = as 

a? 7 



u Pl]\ u lO> aa=l k cr 12 cry crjv-iAr 



(<Ji2)- ■ ■pkl{Ckl)- ■ •PiV-lJv(cjV-lJv) 



+ r In py(o-y) + T + T ^ ^ Pu(<7ki)ln Pki(°ki) + Ay \ = 0. (BIO) 

fei^ij o" fci =l 
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This leads to 



py( ffy )=4e~* E ^~ E V"£-«-"£-*-i^ (Bll) 



Using the normalization condition ( |B6| ), we obtain the factor 



.4 = exp ( - 1 - X] Pn{vki)}& Pu{<Tki) 

kl^ij cr kl = l 



(B12) 



as 



.4 



E 



From Eqs. ( |B11| ) and (|B13| ), the marginal distribution Pijfaj) reads 



Ptji^j) 



(B13) 



y-0 "fE, 12 -Eo-jy -E CTAr _ 1JV W(W)Pl2(o-12)-Pfei(°-fci)-PJV-lJv(o-JV-lJv) 



(B14) 



In order to calculate the sum X)<r 12 ' ' ' E 



(T 12 A^CTjVjV-lV /' 



we rewrite the Hamiltonian 7~t({<j}) 



as 



^({ cr }) = — ^(^ij'+l + ^ij'-l + + — ZhTijCTij + (2 J + h)((Tij] 



(2j + /i) EK) 2 



(B15) 



and using the relations between the local magnetization and the marginal distribution, 
namely, m iJ+1 = YS u+1 =i^ij+ip{^%j+i), etc., we obtain 

-/fS E ^({°'})Pl2(0-12)---Pfe/( 'fe/)---PiV-liv(0-JV-lJv) 

(712 Vkl &N-1N 



J 



2h 



1 



(2 J + h) 
T 



- Tp T~t{{°})pl2{0l2)---Pkl{0kl)---PN-IN{0N- 



1N) 



(B16) 



Me/' 



where i' stands for a set of the sites except for 
Using (|B16| ), we rewrite Pifaij) as 
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Q±[J(mij +1 +mij- 1 +m i+1 j+m i - l! j)+2hTij]<7ij- (2J + ft) (ct^) 2 

Pij((?ij) — -p. Y77r \ \ \ 7777 \ (2j+h) , 77' 

where the factors 

appearing in both numerator and denominator of the Pij{cTij) were canceled. As the 
results, we obtain as follows. 

Q 



m 



'.J 



J2 a ijpij( a i 



(j. . e ^[J(rn i j +1 +m i j- 1 +m i+1 j+m i ^ 1 j)+2hT ij ]a ij -^±^((Ti : j) 2 

^" ' " . (B18) 



'13 



EQ p ^[J(m i j +1 +m i j- 1 +rni +1>j +rni- 1 j)+2hT ij }a i j- ^ 2J r^ h ^ (<Jij) 2 
<Tij=l c 

If we set T = 1, we can obtain the recursion relations with respect to the local magneti- 
zation rriij under the periodic boundary condition as 

J*) (a) 



EKcJ (jj . . . 



= {J[m!$ +1 + mg_! + mg^ + m&j] + a - (2 J + fe)^, (B20) 

m iii+i v = rriij, m+Nj = mj, (B21) 
which were obtained in the previous section as Eqs. (|48D , (|49|) and (|5C|). 
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FIGURES 

FIG. 1. The magnetization of the original image for the case of Q = 3 (a) and Q = 4 (b). 
The solid lines correspond to globally stable solutions. 

FIG. 2. The Hamming distances without glassy term {(3j = 0) for the case of Q = 3 (a)(d). 
The figure (d) is obtained by expanding the figure (a) around its minimum. The magnetization 
m and corresponding free energy — /rs are plotted in (b) and (c), respectively. Here, the dots 
lines and the solid lines are locally stable states and globally stable states, respectively. 

FIG. 3. The Hamming distances without glassy term {(5j = 0) for the case of Q = 4 (a)(d). 
The figure (d) is obtained by expanding the figure (a) around its minimum. The magnetization 
m and corresponding free energy — /rs are plotted in (b) and (c), respectively. 

FIG. 4. The Hamming distance Dr of the MAP estimate for the case of 

Q = 3, T s = 0.75, ro = r = 1.0 (a). L>r is plotted as a function of the scaled field H = h/fa. 
We see that the minimum of L>r is appeared at H = H opt = to/2t 2 (3 s = 0.375. The Hamming 
distance of the MPM estimate is plotted in (b) as a function of the temperature Td for several 
values of H. The figure shows that the minimum of the MPM estimate with H = H opt is lower 
than that of the MAP estimate with H = H opt . 

FIG. 5. The Hamming distance Dr as a function of the strength of the glassy term (5j for 
several values of Jo. We set J = 1.0. The point (3j = corresponds to the minimum of the 
Hamming distance in FIG. 2. 

FIG. 6. The time evolutions of the Hamming distance are plotted in (a) = T s = 0.75 
and (b) Td = 0.2 ^=T S = 0.75. We see that for the case of Td = T s , the Hamming distance 
converges to its optimal value for any initial state of the dynamics. On the other hand, if we 
set Td = 0.2^T S , the Hamming distance converges to the wrong state which is higher than the 
Hamming distance between the original image and the corrupted one, that is, Djp = 0.5. 
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FIG. 7. The inverse relaxation time 1/to as a function of Td for the case of 
Q = 3,T S = 0.75, r = r = 1.0. 

FIG. 8. The Hamming distance calculated by Monte Carlo simulation for 100 x 100 standard 
picture "house". The curves averaged over 20 MCS runs are shown in (j3j = (a) and (3j ^ 
(b))- 

FIG. 9. The original picture (a) ("house", size 100 x 100), the corrupted picture by a = 1.2 
Gaussian noise (b), the restored pictures at (3j = (c) and (3j ^ (d) are displayed 

FIG. 10. The results of the iterative algorithm are displayed. The original "girl" picture of 
128x128 (a), the degraded picture (b), the restored pictures with J = 0.2 (c), J = 1.8 (d) and 
J = 2.5 (e) are shown. 

FIG. 11. The results of the iterative algorithm are displayed. The original "chair" picture of 
128x128 (a), the degraded picture (b), the restored pictures with J = 0.2 (c), J = 1.8 (d) and 
J = 2.5 (e) are shown. 

FIG. 12. The results of the iterative algorithm are displayed. The original "house" picture 
of 128x128 (a), the degraded picture (b), the restored pictures with J = 0.2 (c), J = 1.8 (d) 
and J = 2.5 (e) are shown. 

FIG. 13. The Hamming distance between the original image and the restored one Z?h as 
a function of J [(a)] obtained by the iterative algorithm. The Hamming distance between the 
restored image and the degraded one D$ is shown in (b). 
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